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Abstract 

The Modified Craig-Sneyd (MCS) scheme is a promising splitting 
scheme of the ADI type introduced by In 't Hout & Welfert [Appl. Num. 
Math. 59 (2009)] for multi-dimensional pure diffusion equations having 
mixed spatial-derivative terms. In this paper we investigate the exten- 
sion of the MCS scheme to two-dimensional convection-diffusion equations 
with a mixed derivative. Both necessary and sufficient conditions on the 
parameter 9 of the scheme are derived concerning unconditional stability 
in the von Neumann sense. 
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1 Introduction 



We consider the numerical solution of initial value problems for large systems 
of ordinary differential equations (ODEs), 

U'(t) = F(t,U(t)) (t>0), 17(0) = 0b, (1-1) 

with given vector- valued function F, given initial vector Uq, and unknown vec- 
tors Ufa (for t > 0). Our interest in this paper lies in systems (|l.ip that arise 
from semi-discretization of initial-boundary value problems for two-dimensional 
convection-diffusion equations possessing a mixed spatial-derivative term, 

du 

— = dnu xx + (d 12 + d 21 )u xy + d 2 2U yy + c x u x + c 2 u y . (1.2) 

Here c = fa) and D = (dij) denote a given real vector and a given positive semi- 
definite real matrix, respectively. A main application area of equations of the 
kind (|1.2[) is financial option pricing theory, where mixed derivative terms u xy 
arise naturally since the underlying Brownian motions are usually correlated to 
each other. Extensive details and examples of financial applications are given 
in, for example, the references [9j [101 E] ■ 

For the numerical solution of semi-discrete problems (jl.ip , splitting schemes 
form an effective and popular means, cf. e.g. [SUB]- This paper is devoted to the 
analysis of a recent splitting scheme of the Alternating Direction Implicit (ADI) 
type that has been tailored so as to deal with equations possessing a mixed 
derivative term. Let 9 > be a given fixed parameter. Assume the right-hand 
side function F is decomposed into a sum 

F{t,v) = F (t,v) + F x (t,v) + F 2 {t,v), (1.3) 

where Fq represents the contribution to F stemming from the mixed derivative 
term, and Fj (for j = 1, 2) represents the contribution to F stemming from all 
spatial derivative terms in the j-th spatial direction. Let At > be a given 
time step and define temporal grid points by t n — n ■ At (n = 0, 1, 2, . . .). We 
consider the following splitting scheme for (|l.ip . generating in a one-step fashion 
successive approximations Ui,U 2 , U3, ... to U(ti),U(t 2 ), Ufa), . . . : 
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Method (|1.4[) is called the Modified Craig-Sneyd (MCS) scheme. It has recently 
been introduced, in a slightly more general form, by In 't Hout & Welfert [4]. 
Taylor expansion yields that the MCS scheme has classical order of consistency 
equal to two for any value 9. 

The MCS scheme can be viewed as an extension of the second-order Craig- 
Sneyd (CS) scheme proposed in [T]. The latter scheme, called "iterated scheme" 
in loc. cit., is equivalent to (jl.4[) with parameter value 9 = i. 

A perusal of (|1.4I) shows that the Fo term is always treated explicitly, whereas 
the Fi and F 2 terms are treated implicitly. More precisely, the MCS scheme 
starts with an explicit Euler step applied to the full system (11.11) which is suc- 
ceeded by two implicit corrections corresponding to each of the two spatial 
directions. Subsequently, an explicit update is performed, which is followed 
again by two implicit, unidirectional corrections. Accordingly, the MCS scheme 
retains the well-known key advantage of ADI schemes over standard implicit 
methods, such as the Crank-Nicolson scheme, that the (linear or nonlinear) 
systems to be solved in each time step are much easier to handle. 

The adaptation of ADI schemes to convection-diffusion equations with mixed 
derivative terms has been studied by a number of authors. Several stability 
results, in the sense of von Neumann, have been obtained. McKcc ct al. [6, 
[7] considered a simpler version of (|1.4I) . which is equivalent to the first two 
lines with U n = Y 2 . This basic scheme, also known as the Douglas scheme, 
is of order one for any value 9 in the presence of a mixed derivative term. 
McKee et al. showed that if 9 = i, then it is unconditionally stable when 
applied to a standard finite difference (FD) discretization of (|1.2[) . Next, Craig 
& Sneyd [T] formulated the second-order CS scheme and proved that this scheme 
is unconditionally stable in the case of (|1.2I) with c = 0. Recently In 't Hout & 
Welfert [3J S] extended the above stability results in various ways. We state here 
the main results pertinent to the situation at hand. Firstly, for the CS scheme 
unconditional stability was proved [3] in the general case of (|1.2p . Secondly, it 
was shown [3] that in the case of (ll.2j) with c = the MCS scheme (|1.4p is 
unconditionally stable whenever 9 > |. 

Up to now it is an important open question when the MCS scheme, with 
9 =^ | , is unconditionally stable in the application to general equations (|1.2[) , 
i.e., with arbitrary c and positive semi-definite D. As it turns out, an analysis 
of this is not straightforward, related to the fact that the eigenvalues of the 
semi-discrete linear operators move from the real line in the pure diffusion case 
to the complex plane in the general, convection-diffusion case. In the present 
paper we shall arrive at positive results on the above question. 

For the stability analysis we consider the linear scalar test equation 



with complex constants Xj (0 < j < 2). When applied to f| 1 . 5[) . the MCS scheme 
(II. 4[) reduces to the scalar iteration 



U'(t) = (Ac + Ai + X 2 )U(t) 



(1.5) 



U n = Sg{z ,Z 1 ,Z 2 ) U n -\ 



(1.6) 
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with z 3 = At ■ Xj (0<j< 2) and 



S e (zo 7 z 1 ,z 2 ) = 1 + 



Z Q + z 



V 



+ 








p2 



(1.7) 



where we use the notation 



z = z\ + Z2 and p = (1 — 0zi)(l 



The iteration (11.61) is stable if 



\Se(z , Zi, z 2 )\ < 1. 



(1.8) 



In the von Neumann framework, the Xj represent eigenvalues of the linear 
operators Fj that are obtained after semi-discretization, on a uniform spatial 
grid, of the convection-diffusion equation (|1.2j) with constant coefficients and 
periodic boundary condition. Corresponding to the positive semi-dehniteness of 
the diffusion matrix D, it was shown in [3] (cf. also Sect. \S§ that for standard FD 
discretizations the following condition on the scaled eigenvalues Zj is fulfilled, 



where all bounds are sharp. In view of this, a natural stability requirement on 
the scheme (|1.4|) when applied to equations (|1.2p with mixed derivative terms 
is that (ll.8[) holds whenever (I1.9P is satisfied. 

An outline of the rest of this paper is as follows. In Sect. [2] we study for 
which parameter values 6 the implication (jl.9[) => (ll.8[) is fulfilled. Four cases 
are investigated, depending on whether zo is real or complex valued and whether 
Z\, Z2 are (both) real or complex valued. In Sect. [3] the results of Sect. [2] are 
applied and discussed relevant to an application of the MCS scheme (ll.4[) to 



2 Stability results for the MCS scheme 

Let / denote the imaginary unit. In this section we study the stability require- 
ment (|1.9p (|1.8I) . The following introductory result gives a criterion on 8 for 
the case Zq — 0. This is pertinent to the situation where no mixed derivative 
term is present in (|1.2j) . 

Theorem 2.1 There holds \Sg(0, Zi, z 2 )\ < 1 for all Z\,zi € C with Sftzi < 0, 
5Rz2 < if and only if 9 > 4 . 

Proof The rational function Sg(0, z\, z 2 ) has no poles in the set 3?zi,5Rz2 < 
and therefore attains its maximum on the boundary of this set. Thus assume 
z\ = Ib\, z 2 = Ib 2 with 61,62 G R- We have 



\z Q \ < 2^/Uz^Mz^, 5Rzi<0, 5Rz 2 <0, 



(1.9) 
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Write u = 1 - 8 2 b 1 b 2 and v = h + b 2 . Then 

p = u — 9Iv , p 2 = u 2 — 8 2 v 2 — 29Iuv , pz — 8v 2 + luv , z 2 = —v 2 
and it follows after some algebraic manipulations that 

\ P 2 +pz + (±- 6)z 2 \ 2 - \p 2 \ 2 = [{6 2 -29 + i) 2 - 

Hence, 

|5e(0, zi, 22)| < 1 whenever ^Rzi,3tz 2 < 

if and only if 



v 4 . 



2* + §1 



which is equivalent to 9 > \. 



In [1] the stability of ADI schemes for pure diffusion equations with mixed 
derivatives was analyzed. This concerns the case where all Zj are real-valued. 
For the MCS scheme and two spatial dimensions, the following criterion on 9 
was obtained. 

Theorem 2.2 There holds \Se(zo, z\, z 2 )\ < 1 whenever zq,zi,z 2 G K satisfy 
(TOj) if and only if9>\. 

Proof See d Thm. 2.5]. ■ 

In most applications, also a convection term is present. Accordingly, one is 
led to considering complex- valued z\,z 2 . The next theorem gives a necessary 
condition on 9 for this situation. 

Theorem 2.3 Suppose \Sg(zo, zi, z 2 )\ < 1 for all zq € K and zi,z 2 £ C satis- 
fying ((213). Then 9> f. 



Proof The result is obtained by a Taylor expansion at the point zo = z\ = 
Z2 = 0. We take zo = —2a and z\ = z 2 — an where n = I + 1 and a G R with 
a j" 0. This choice was found to be convenient after numerical experimentation. 
Inserting into HUT]) and using 1/(1 -£) = l+ £ + £ 2 + (f -> 0), it follows 

that 

5 e (-2a, af7 , a?7 ) = l + _^ + (tf-I-^)^-^ 

= 1 + 21 a - 2a 2 + (209 2 - 89 - 89I)a 3 + C(a 4 ). 

This yields 

|S' (? (-2a,a?7,ar/)| 2 = 1 + (406» 2 - I69)a 3 + C(a 4 ). 

The right-hand side is bounded by 1 for a f only if A09 2 — 168 > 0. Hence, it 

must hold that 6 > |. ■ 

— 5 
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Based on strong numerical evidence (see Sect. [3]) we conjecture that the condi- 
tion on 9 in Theorem 12. 31 is also sufficient, but a proof is currently lacking. 

The above results dealt with real- valued z . The two subsequent theorems 
concern arbitrary, complex- valued zq. A preliminary result is 

Lemma 2.4 Let a, b,c G R be given. If \a + b + c\ = l and \a( 2 + b( + c\ < 1 
whenever £ G C wrei/i |£| = 1, i/ien a6 + 6c + 4ac > 0. 

Proof Consider the function / defined by 

f(ip) = \ae 2Iv + be Iv + c\ 2 ((p G R). 

There holds 

/(<p) = [a cos(2^) + 6 cos(^) + c] 2 + [a sin(2</j) + b sin(^)] 2 
= a 2 + b 2 + c 2 + 2abcos(ip) + 2bccos(tp) + 2accos(2<y9). 

One readily verifies that /(0) = 1, /'(0) = 0, /"(0) = -2(ab + be + 4ac) and 
hence 

f(tp) = l-{ab + bc + 4ac)tp 2 + 0(tp 3 ) (<p -> 0). 
Using that f(ip) < 1 whenever ip G R, proves the assertion. ■ 

For the case where Zi, Z2 are real- valued, we obtain the following necessary lower 
bound on 9. Numerical experiments indicate that this bound is sufficient as well. 

Theorem 2.5 Suppose \S$(zq, z\, z%)\ < 1 for all zq G C and z\,z% G R satis- 
fying fO)l . JTiera 0> ^. 

Proof Setting g = p 2 + + — #)z 2 and to = p + (1 — 9)z, we can write 

c i \ \ z l + wz ° + g fo ml 

Let y = 2y / zi^2- Since |iSg(,zo, zi, 22)! < 1 f° r ah z o G C with \zq\ < y we have 

|iy 2 C 2 + wyC + «| < P 2 for all C G C with |C| < 1. 

Assume z\ = zi. Then z = —y and it is easily seen that \y 2 + wy + q = p 2 . 
Therefore Lemma l2T4"l can be applied and, using y > 0, this leads to the necessary 
condition 

\wy 2 + qw + 2qy > whenever Z\ = z-i < 0. (2-H) 
Denote x = 9y. Then 

p = 1 + x + \x 2 . 
Next, after some computations, there follows 

\wy 2 +qw + 2qy=p 3 +p 2 x - hx 3 + 2px 2 ). 



G 



By (|2.1ip . we arrive at 

x 3 + 2px 2 

V > — 5 o — • 

pC _|_ pZj, 

The right-hand side is a rational function of x > 0, which is readily seen to have 
a global maximum at x — 2. Inserting this value yields the lower bound 9 > yrj. 



The final result in this section concerns the most general case, where all Zj 
are complex- valued. To derive this result we employ a lemma from [3] pertinent 
to the condition (|1.9j) . For completeness, its concise proof is included here. 



Lemma 2.6 7/z 1 ,z 2 G C with 3tzi < 0, $lz 2 < 0, then 



p 






26 




29 +Z 



Proof Define the vectors 



Their Euclidean norms are 



yJ-Tkz.j 
\l + 6zj\/V2d 



J = 1,2. 



^Il = \ -2^ + 



29 



Next, their standard inner product is 
(vi, v 2 

Applying the Cauchy-Schwarz inequality gives 



2V^ + " 1 + fc ''' 1+fc '" 
29 



which concludes the proof. 



P_ 

29 



< 



29 

2y/Wz^Wz~ 2 
-9z 2 \ 



29 



+ z 



29 



For the most general case, we have the following positive result: 

Theorem 2.7 If | < < 1, then \Sg(z , z\, z 2 )\ < 1 whenever z ,z\,z 2 € C 
satisfy (fOj) . 



Proof The expression (12.10)) for So yields 

2 



\Se(z , zi,z 2 )\ < i 



z 

1 + - 



By invoking Lemma 12.61 it follows that \Sg(zo, z\, z 2 )\ is bounded from above by 



2\29 



1 

29 



1 

29 



1 

20 



i + (i-ey- 



z 
V 
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We can write 

1 + 29- = re Iv with < r < 1 and < tp < 2ix. 
P 

Define 

h(tp,r) = \26 + (l-6)(re^-l)\, 

f 2 (<p,r) = \W 2 + 49(re Iip - 1) + (1 - 26)(re Iip - 1) 2 | . 

Then it follows that 

Let | < 9 < 1. We prove that the right-hand side of (|2.12[) is bounded 
by 1 for all < r < 1, < tp < 2n. First note that fj{2n — ip,r) = fj(<p,r) 
(j = 1, 2) and therefore it suffices to consider < ip < tt. Let r G [0, 1] be fixed 
but arbitrary and define cjj(<p) — fj((p,r) 2 (j = 1,2). For the function gi it is 
readily verified that 

g[(ip) = -2(36»- 1)(1 - 6)r sin^. 

This directly implies that gi, and hence fi, is nonincreasing on [0, tt]. For the 
function gi a more elaborate computation shows 

g' 2 ((p) = 4(29 - 1) (4:9 - 1) [2(26* - l)r cos tp + r 2 - (49 - 1)] rsin^s. 

In view of 

2(20 - l)r cos tp + r 2 - (49 - 1) < 2(29 - 1) + 1 - (49 - 1) = 

we find that also g 2 , and hence f 2 , is nonincreasing on [0,7r]. Consequently, it is 
sufficient to prove that the right-hand side of (|2.12j) is bounded by 1 whenever 
< r < 1, <p = 0. Write s — r — 1 e [—1, 0]. One easily verifies that 

A(0,r) = 29 + {1-9)8, 

/ 2 (0,r) = 86 ,2 + 46» s + (l-26 , )s 2 . 

Inserting this and rearranging terms, it follows that the upper bound (|2.12[) is 
(in fact) equal to 1 whenever < r < 1, (p = 0. m 

Numerical evidence leads to the conjecture that the conclusion of Theorem 12.71 
is valid for all 9 > j^, i.e., under the (necessary) lower bound of Theorem 12.51 
A proof of this does not appear to be straightforward. We note that in the proof 
above the assumption i < 9 < 1 is used in an essential manner. 
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3 Application and discussion 



In this section we discuss an application to convection-diffusion equations (|1.2[) . 
We semi-discretize on the unit square [0, 1] x [0, 1] by using central second-order 
FD schemes on a Cartesian grid with mesh widths Ax and Ay in the x and y 
directions, respectively: 

W U « " " 2A ; " (3-la) 

- " " 2 A//' ' (3 - lb) 



- ^ (A^ (3 - lc) 
(««/)*.,• - (3-ld) 



(u X y) 



VJi,j 



(1 + + Mj-lj-i) ~ (1 ~ gX^-lJ+l + Ui+lJ-l) 

4AxAy 



(3.1e) 



4AxAy 

Here /3 denotes a real parameter with — 1 < j3 < 1 and we use the notation 
Uij = u(iAx,jAy,t). We note that the right-hand side of (|3.1e ) is the most 
general form of a second-order FD approximation of the mixed derivative u xy 
based on a centered 9-point stencil. When (3 = 0, it reduces to the well-known 
4-point formula 

i x + Uj-l.j-l — "i-lj'+l — Ui+l,j-l 

[Uxy) ^ ~ TaTa^ ' 

Assuming constant coefficients and a periodic boundary condition for (|1.2[) , 
the above FD discretization yields a splittcd, semi-discrete system (jl.ip . (II .3[) 
where Fj(t,v) = Aj*u for j = 0,1,2 with constant matrices A, . The matrix Aq 
represents the cross derivative term in (|1.2[) and At, A2 represent the spatial 
derivatives in the x and y directions, respectively. The periodicity condition 
implies that the Aj are Kronecker products of circulant (thus normal) matrices 
that commute with each other, and are therefore simultaneously diagonalizable 
by a unitary matrix. Hence, stability can be rigorously analyzed by considering 
the scalar test equation (|1.5p with A., eigenvalues of Aj (0 < j < 2). This is 
equivalent to a von Neumann stability analysis. By inserting discrete Fourier 
modes, it follows that the scaled eigenvalues Zj are given by 



zi = 

Z2 

where 



— — — — . — ...j 

= (d 12 + d 2 \) b [- sin 0i sin0 2 + - cos0i)(l - cos0 2 )] , (3.2a) 

= — 2dnai(l — cos 0i ) + Ic\qi sin 0i , (3.2b) 

= -2d 22 a 2 (l - cos0 2 ) + Ic 2 q 2 sin 02 , (3.2c) 



At At At _ At _ At 

ai_ (A^' a2 "(A^' ~~ ~AxAy ' qi ~ ~Ax ' 92 ~~ Ay 
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The angles <pj are integer multiples of 27r/m 3 (j = 1,2) where mi, mi are the 
dimensions of the grid in the x and y directions, respectively. 

Using the positive semi-definiteness of the diffusion matrix D, an elementary 
calculation shows [3] that zo,zi,Z2 fulfill the condition (jl.9[> . independently of 
At, Ax, Ay. Upon invoking Theorem 12.71 the following neat stability result is 
obtained for the MCS scheme applied to (jl.2p . 

Theorem 3.1 Consider equation with positive semi-definite matrix D and 
periodic boundary condition. Let the semi-discrete system U.l\) . \1.3\) be ob- 
tained after FD discretization and splitting as described in this section. Then 
the MCS scheme \l-4ty is unconditionally stable when applied to U.l\) . U.3\) 
whenever i < 6 < 1. Moreover, this conclusion remains valid when any other 
stable FD discretizations for u x , u y are used in place of HS.lh ). 13.1b ). 

The last part of Theorem 13.11 follows directly from the fact that the real parts 
of the new eigenvalues z\, z^ are always smaller than those of (|3.2b ). (|3.2b ). 
respectively, and hence, (|1.9I) remains true. 

An inspection of p. 2a ) yields that the eigenvalues zq have the property that 
their imaginary part is identically equal to zero. Accordingly, it is of particular 
interest to know all parameter values 9 such that the stability requirement 
(|1.9P =>■ (|1.8[) holds for just real- valued zq. Theorem 12.71 provides the sufficient 
condition i < 6 < 1, whereas Theorem [2?3] yields the necessary condition 9 > |. 

Next, we remark that the MCS scheme has recently been applied successfully 
in [2] to actual convection-diffusion equations (|1.2p with mixed derivative terms 
using the parameter value 9 = g. This seems to be surprising, as this value 
was determined |H for pure diffusion equations (|1.2I) and it clearly does not 
satisfy the necessary condition 9 > | for equations with convection. We note 
that reasons for choosing a smaller 9 in the MCS scheme are a reduced error 
constant and better damping properties compared to the original CS scheme, 
see [2]. 

Theoretical results on the latter two issues are not known at this moment. 
To gain insight, we have performed a numerical experiment. Let rj. and r^j for 
i,j = 1,2 denote independent, uniformly distributed random numbers in [0, 1] 
and consider random triplets (zq, z\,z-i) given by 

z = (2r M - 1) • 2y/KziKz2 and Zj = -lO 1 " 5 ^ ± / lO 1 " 5 ^--* (7 = 1,2). 

Then (TTU) holds and z € M. For each 9 = \ + ^ with k = 0, 1, . . . , 100 we com- 
puted the maximum value of \Se(zQ, z\, z^)\ over two million points (zo, Z\, Z2) 
above. The outcome is displayed in Figure Q] 

Figure [1] reveals the intriguing result that the estimated maximum value of 
I Sg I is very close to 1 whenever 9 > | . For 9 = | we arrive at a maximum value 
of 1.02. Additional experiments in this case suggest that \Sg \ is larger than 1 for 
a limited set of points (zq, z\,Z2), and at most 1 under only a slightly stronger 
condition on zq than in (11.91) . Because of these observations, it is very plausible 
that the MCS scheme performs well in actual applications to (|1.2p . also with 
convection, already when 9 = ^. 
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Figure 1: Estimated maximum of \Sg(zo, z\, Z2)\ under (|1.9[) with zq G 



Subsequently, an examination of the obtained numerical results indicates 
that \Sg\ < 1 for all 9 > %. This supports our conjecture formulated below 
Theorem [231 

In view of the above, it is likely that the condition on 9 in Theorem 13. II can 
be relaxed to 9 > |, and next, that a slightly modified version of Theorem 13. II 
holds under the (weaker) assumption 9 > |. In future research we intend to 
study these issues theoretically. 
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